function [tab_BS] = pool_BS(s_now_batch,rule,varargin) 
% This function #pools# bootstrap results from a designated folder with a
% specific file name prefix

% Input: 
% (1) s_now_batch: time stamp of bootstrap initilization 
% (2) rule: quadrant or cubic, different fields are extracted based on rule
% (3) (optional) fpath_BS: path to folder where individual boot.mat is
%                           saved
% Output 
% (1) tab_BS: table of bootstrap results (each row in the table represent a
% specific bootstrap run)

%% parse varargin
if nargin==2
    fpath_BS = [];
else
    fpath_BS = varargin{1};
end
    
%% Specify common filename prefix
s_fformat = sprintf('%s*_BS_*.mat',s_now_batch);
%% Choose FIles
if isempty(fpath_BS)
    % Manually select files to load ---------
    [c_fbase_BS,fpath_BS] = uigetfile(s_fformat,'multiselect','on');
    if ~iscell(c_fbase_BS)
        c_fbase_BS={c_fbase_BS};
    end
else
    % Automatically select all files in the folder
    oldpath=pwd;
    cd(fpath_BS);
    
    a=dir(s_fformat);
    c_fbase_BS = {a.name}';
    cd(oldpath);
end
N_file = length(c_fbase_BS);
%% Verify BS iterations
expression='(?<=BS_)\d+(?=.mat)';
aBSidx=cellfun(@(x) str2double(regexp(x,expression,'match','once')),c_fbase_BS','uni',1);

%-- Are there any duplicative values
[uBSidx,~,ic]=unique(aBSidx,'stable');
bit_dup=zeros(size(uBSidx));
for i=1:length(uBSidx)
   bit_dup(i)=(sum(ic==i)>1);
end
if any(bit_dup)
    fprintf('Duplicated BS idx :{%s}',num2str(uBSidx(bit_dup)))
end

%% Load
struct_BS=struct();
for ifile = 1:N_file
    fprintf('pooling No.%3.0f/%3.0f\n',ifile,N_file);
    fbase = c_fbase_BS{ifile};
    fpathf = fullfile(fpath_BS,fbase);
    if rule=="quadrant" || rule=="onedim"
        % --- Load
        S=load(fpathf,'bsi','s_now_batch','Wb','maxWb','minWb',...
            'scaled_ATT_bs','ATE_bs','W_bs_Ghat','bsperm','dt','covariate');
        % --- organize
        struct_BS(ifile).bsi=S.bsi;
        struct_BS(ifile).Wb = S.Wb;
        struct_BS(ifile).maxWb = S.maxWb;
        struct_BS(ifile).minWb = S.minWb;
        struct_BS(ifile).scaled_ATT_bs = S.scaled_ATT_bs;
        struct_BS(ifile).ATE_bs = S.ATE_bs;
        struct_BS(ifile).W_bs_Ghat = S.W_bs_Ghat;
        struct_BS(ifile).bsperm = S.bsperm;
        struct_BS(ifile).dt = S.dt;
    elseif rule=="cubic"
        % --- Load
        S=load(fpathf,'bsi','s_now','v','v_pd','v_pi','v_nd','v_ni','in_Ghat','g_hat',...
            'scaled_ATT_bs','ATE_bs','W_bs_Ghat','bsperm','dt');
        % --- organize
        struct_BS(ifile).bsi=S.bsi;
        struct_BS(ifile).v = S.v;
        struct_BS(ifile).v_pd = S.v_pd;
        struct_BS(ifile).v_pi = S.v_pi;
        struct_BS(ifile).v_nd = S.v_nd;
        struct_BS(ifile).v_ni = S.v_ni;
        struct_BS(ifile).in_Ghat = S.in_Ghat;
        struct_BS(ifile).g_hat = S.g_hat;
        struct_BS(ifile).scaled_ATT_bs = S.scaled_ATT_bs;
        struct_BS(ifile).ATE_bs = S.ATE_bs;
        struct_BS(ifile).W_bs_Ghat = S.W_bs_Ghat;
        struct_BS(ifile).s_now_shared = S.s_now;
        struct_BS(ifile).bsperm = S.bsperm;
        struct_BS(ifile).dt = S.dt;
    end
end
tab_BS = struct_BS;

end